###Replication File for: Blake, Daniel J. 2013. Thinking Ahead: Government Time Horizons and the Legalization of International Investment Agreements

library(foreign)
library(boot)
library(MASS)
library(logistf)
library(sandwich)
library(lmtest)
library(Design)
library(pscl)

setwd ("C:\\Users\\Daniel Blake\\Dropbox\\IO FIles\\Replication")
#setwd("/Users/danielblake/Dropbox/IO FIles/Replication")

#Read in the data

authost <- read.dta("carveouts_hostaut.dta")
demhost <- read.dta("carveouts_hostdem.dta")
authostdem <- read.dta("carveouts_hostauthomedem.dta")
demhostdem <- read.dta("carveouts_hostdemhomedem.dta")

#Attach the Data

#Autocratic Host
attach(authost)

authost.nbreg<-as.data.frame(cbind(ntcarveouts, host_surv , nthsth 
	, host_yearsindep , wealthratio , yearcount,  host_xconst 
	, host_gdppc_growth , hostcode_2006, homecode_2006,  home_xconst 
	, home_dem))
authost.nbreg<-na.omit(authost.nbreg)
dim(authost.nbreg)
class(authost.nbreg)

authost.sel<-as.data.frame(cbind(ntcarveouts, host_surv , nthsth 
	, host_yearsindep , wealthratio , yearcount,   home_xconst , host_xconst  
	, host_gdppc_growth , hostcode_2006, homecode_2006,  pr_host_c1, pr_home_c1, home_dem))
authost.sel<-na.omit(authost.sel)
dim(authost.sel)
class(authost.sel)

authost.hurdle<-as.data.frame(cbind(dv_hurdle,  host_surv , nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst , home_xconst  
	, host_gdppc_growth , hostcode_2006, homecode_2006, home_dem))
authost.hurdle<-na.omit(authost.hurdle)
dim(authost.hurdle)
class(authost.hurdle)

detach(authost)

##Autocratic Host, Democratic Home

attach(authostdem)
authostdem.nbreg<-as.data.frame(cbind(ntcarveouts, host_surv, home_govagelp, home_agedem, nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst 
	, host_gdppc_growth , hostcode_2006, homecode_2006,  home_xconst ))
authostdem.nbreg<-na.omit(authostdem.nbreg)
dim(authostdem.nbreg)
class(authostdem.nbreg)

authostdem.sel<-as.data.frame(cbind(ntcarveouts, host_surv, home_govagelp, home_agedem, nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst , home_xconst   
	, host_gdppc_growth , hostcode_2006, homecode_2006,  pr_host_c2, 	pr_home_c2))
authostdem.sel<-na.omit(authostdem.sel)
dim(authostdem.sel)
class(authostdem.sel)


authostdem.hurdle<-as.data.frame(cbind(dv_hurdle, host_surv, home_govagelp, home_agedem, nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst , home_xconst  
	, host_gdppc_growth , hostcode_2006, homecode_2006))
authostdem.hurdle<-na.omit(authostdem.hurdle)
dim(authostdem.hurdle)
class(authostdem.hurdle)

detach(authostdem)

#Democratic Host

attach(demhost)

demhost.nbreg<-as.data.frame(cbind(ntcarveouts, host_govagelp, host_agedem,  nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst 
	, host_gdppc_growth , hostcode_2006, homecode_2006,  home_xconst , home_dem))
demhost.nbreg<-na.omit(demhost.nbreg)
dim(demhost.nbreg)
class(demhost.nbreg)

demhost.sel<-as.data.frame(cbind(ntcarveouts, host_govagelp, host_agedem,  nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst , home_xconst   
	, host_gdppc_growth , hostcode_2006, homecode_2006,  pr_host_c3, 	pr_home_c3, home_dem))
demhost.sel<-na.omit(demhost.sel)
dim(demhost.sel)
class(demhost.sel)

detach(demhost)

#Democratic Host, Democratic Home

attach (demhostdem)

demhostdem.nbreg<-as.data.frame(cbind(ntcarveouts, host_govagelp, host_agedem,  nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst 
	, host_gdppc_growth , hostcode_2006, homecode_2006, home_govagelp, home_agedem,   home_xconst ))
demhostdem.nbreg<-na.omit(demhostdem.nbreg)
dim(demhostdem.nbreg)
class(demhostdem.nbreg)

demhostdem.sel<-as.data.frame(cbind(ntcarveouts,  host_govagelp, host_agedem,  nthsth 
	, host_yearsindep , wealthratio , yearcount,   host_xconst , home_xconst   
	, host_gdppc_growth , hostcode_2006, homecode_2006, home_govagelp, home_agedem,   pr_host_c4, 	pr_home_c4))
demhostdem.sel<-na.omit(demhostdem.sel)
dim(demhostdem.sel)
class(demhostdem.sel)

detach(demhostdem)


#Functions for clustering standard errors 
#Adapted from code written by Mahmood Arai (Stockholm University) and presented in
#his research note "Cluster-robust standard errors using R"

 cl   <- function(dat,nbreg, cluster){
              attach(dat, warn.conflicts = F)
              library(sandwich)
              M <- length(unique(cluster))   
              N <- length(cluster)  	        
              K <- nbreg$rank		             
              dfc <- (M/(M-1))*((N-1)/(N-K))  
              uj  <- apply(estfun(nbreg),2, function(x) tapply(x, cluster, sum));
              vcovCL <- dfc*sandwich(nbreg, meat=crossprod(uj)/N)
              coeftest(nbreg, vcovCL) }


  mcl <- function(dat,nbreg, cluster1, cluster2){
              attach(dat, warn.conflicts = F)
              library(sandwich);library(lmtest)
              cluster12 = paste(cluster1,cluster2, sep="")
              M1  <- length(unique(cluster1))
              M2  <- length(unique(cluster2))   
              M12 <- length(unique(cluster12))
              N   <- length(cluster1)          
              K   <- nbreg$rank             
              dfc1  <- (M1/(M1-1))*((N-1)/(N-K))  
              dfc2  <- (M2/(M2-1))*((N-1)/(N-K))  
              dfc12 <- (M12/(M12-1))*((N-1)/(N-K))  
              u1j   <- apply(estfun(nbreg), 2, function(x) tapply(x, cluster1,  sum)) 
              u2j   <- apply(estfun(nbreg), 2, function(x) tapply(x, cluster2,  sum)) 
              u12j  <- apply(estfun(nbreg), 2, function(x) tapply(x, cluster12, sum)) 
              vc1   <-  dfc1*sandwich(nbreg, meat=crossprod(u1j)/N )
              vc2   <-  dfc2*sandwich(nbreg, meat=crossprod(u2j)/N )
              vc12  <- dfc12*sandwich(nbreg, meat=crossprod(u12j)/N)
              vcovMCL <- vc1 + vc2 - vc12
              coeftest(nbreg, vcovMCL)}

##########################################################
##Table 1 -Negative Bionmial regression of NT Carve Outs##
##########################################################


#Table 1 Column 1

nbreg<- glm.nb(ntcarveouts ~ host_surv + home_dem + nthsth 
	+ host_yearsindep + wealthratio + host_gdppc_growth + yearcount 
	+ host_xconst + home_xconst , data=authost.nbreg)

#summary(nbreg)

mcl(authost.nbreg,nbreg, hostcode_2006, homecode_2006) # Two-way clustering on host and home state 			

#cl(authost.nbreg,nbreg, homecode_2006) 		   # Clustering on home state 


#Table 1 Column 2

nbreg <- glm.nb(ntcarveouts ~ host_surv + home_govagelp 
	+ home_agedem + nthsth 
	+ host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst ,
	 data=authostdem.nbreg)

#summary(nbreg)

mcl(authostdem.nbreg,nbreg, hostcode_2006, homecode_2006) # Two-way clustering on host and home state 	

#cl(authostdem.nbreg,nbreg, homecode_2006)   # Clustering on home state 


#Table 1 Column 3

nbreg <- glm.nb(ntcarveouts ~ host_govagelp + host_agedem 	+ home_dem + nthsth + host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst , data=demhost.nbreg)

#summary(nbreg)

mcl(demhost.nbreg,nbreg, hostcode_2006, homecode_2006) # Two-way clustering on host and home state

#cl(demhost.nbreg,nbreg, homecode_2006)	  # Clustering on home state

#Table 1 Column 4

nbreg <- glm.nb(ntcarveouts ~ host_govagelp + host_agedem
	+ home_govagelp + home_agedem + nthsth 
	+ host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst ,
	 data=demhostdem.nbreg)

#summary(nbreg)

mcl(demhostdem.nbreg,nbreg, hostcode_2006, homecode_2006) # Two-way clustering on host and home state

#cl(demhostdem.nbreg,nbreg, homecode_2006)	# Clustering on home state


##########################################################
#########Table 2 - Hurdle Models of NT Carve Outs#########
##########################################################


#Table 2 Columns 1&2

hurdle<- hurdle(dv_hurdle ~ host_surv + home_dem + nthsth 
	+ host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst , data=authost.hurdle, dist = "negbin")

summary(hurdle)

#Table 2 Columns 3&4

hurdle <- hurdle(dv_hurdle ~ host_surv + home_govagelp 
	+ home_agedem + nthsth 
	+ host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst ,
	 data=authostdem.hurdle, dist = "negbin")

summary(hurdle)

##########################################################
###Table 3 -Negative Bionmial Controlling for Selection###
##########################################################


#Table 3 Column 1

authost.nbreg.boot <- function(data, indices){
             
             data<-data[indices,]
             testt <- glm.nb(ntcarveouts ~ host_surv + nthsth
	+ host_yearsindep + wealthratio + host_gdppc_growth 
	+ yearcount + host_xconst + home_xconst + pr_host_c1 + pr_home_c1, data=data)
             #Returns coefficients vector 
             coefficients(testt)
             }


#Return coefficient estimates with boostratpped standard errors

bootstrap<-boot(data=authost.sel, statistic=authost.nbreg.boot, 2000)             

bootstrap

##Generate bias corrected and accelerated confidence intervals for estimates of statistical significance

#90 percent confidence intervals

boot.ci(bootstrap, conf = 0.90, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=11, type=c("bca"))

#95 percent confidence intervals

boot.ci(bootstrap, conf = 0.95, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=11, type=c("bca"))

#Table 3 Column 2

authostdem.nbreg.boot <- function(data, indices){
             
             data<-data[indices,]
             testt <- glm.nb(ntcarveouts ~ host_surv + home_govagelp 
             + home_agedem + nthsth + host_yearsindep + wealthratio 
             + host_gdppc_growth + yearcount + host_xconst + home_xconst  + pr_host_c2 + pr_home_c2, data=data)
             #Returns coefficients vector 
             coefficients(testt)
             }
       
#Return coefficient estimates with boostratpped standard errors
                         
bootstrap<-boot(data=authostdem.sel, statistic=authostdem.nbreg.boot, 2000)  

bootstrap

##Generate bias corrected and accelerated confidence intervals for estimates of statistical significance

#90 percent confidence intervals

boot.ci(bootstrap, conf = 0.90, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=12, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=13, type=c("bca"))

#95 percent confidence intervals

boot.ci(bootstrap, conf = 0.95, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=12, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=13, type=c("bca"))

#Table 3 Column 3

demhost.nbreg.boot <- function(data, indices){
             
             data<-data[indices,]
             testt <- glm.nb(ntcarveouts ~ host_govagelp  
             + host_agedem + nthsth + host_yearsindep 
             + wealthratio 
             + host_gdppc_growth + yearcount + host_xconst + home_xconst  + pr_host_c3 
             + pr_home_c3, data=data)
             #Returns coefficients vector 
             coefficients(testt)
             }
             
#Return coefficient estimates with boostratpped standard errors
             
bootstrap<-boot(data=demhost.sel, statistic=demhost.nbreg.boot, 2000)             

bootstrap

##Generate bias corrected and accelerated confidence intervals for estimates of statistical significance

#90 percent confidence intervals

boot.ci(bootstrap, conf = 0.90, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=12, type=c("bca"))

#95 percent confidence intervals

boot.ci(bootstrap, conf = 0.95, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=12, type=c("bca"))

#Table 3 Column 4

demhostdem.nbreg.boot <- function(data, indices){
             
             data<-data[indices,]
             testt <- glm.nb(ntcarveouts ~ host_govagelp + host_agedem + home_govagelp 			 + home_agedem + nthsth + host_yearsindep + wealthratio 
             + host_gdppc_growth + yearcount + host_xconst + home_xconst  + pr_host_c4 
             + pr_home_c4, data=data)
             #Returns coefficients vector 
             coefficients(testt)
             }

#Return coefficient estimates with boostratpped standard errors
             
bootstrap <- boot(data=demhostdem.sel, statistic=demhostdem.nbreg.boot, 2000)             

bootstrap

##Generate bias corrected and accelerated confidence intervals for estimates of statistical significance

#90 percent confidence intervals

boot.ci(bootstrap, conf = 0.90, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=12, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=13, type=c("bca"))
boot.ci(bootstrap, conf = 0.90, index=14, type=c("bca"))

#95 percent confidence intervals

boot.ci(bootstrap, conf = 0.95, index=1, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=2, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=3, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=4, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=5, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=6, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=7, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=8, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=9, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=10, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=11, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=12, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=13, type=c("bca"))
boot.ci(bootstrap, conf = 0.95, index=14, type=c("bca"))
